AD-A062  165  NAVAL  POSTGRADUATE  SCHOOL  MONTEREY  CALIF  F/G  12/1 

A SIMULATION  STUDY  OF  A CLASS  OF  FIRST-COME  FIRST-SERVED  QUEUES— ETC (U) 
SEP  78  P BOONSONG 


i iun  accTcrcn  wi 


S9IS9  0VW  woo  31IJ  00G 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


D D Ov 


Mfiifk 


DEC  15  19(78 


iLiiawnymU 


HESIS 


r 


P.A.W.  Lewis 


Approved  for  public  release;  distribution  unlimited. 

78  12  11  I9i> 


SECURITY  CLASSIFICATION  OP  THIS  PASS  nw>—  Oats  SlNnO 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


ATALOO  NUMOKN 


«.  titlc  r«f  jhammaj 

A Simulation  Study  of  a Class 
Come  First-Served  Queues  with 
Correlation  Structure 

of  First- 
EARMA 

J.  AUTHOR?  •> 

Prasert  Boonsong 

/ 

S.  TYPE  OP  REPORT  A PERIOD  COVERED 

Master's  Thesis; 
September  1978 


• I.  CON  TOOL  LINO  OPPICS  NANS  ANO  AOONKSS 

Naval  Postgraduate  School 
Monterey,  California  93940 


NCY  NANS  A AC 


SSflf  «Mn*l  i 


i CmmI fin.  Otti*9j 


It.  REPORT  OATS 

September  1978 


IS.  SKCURITY  CLASS.  (•!  till  a 

Unclassified 


ft.  OtSTftlEUTlOM  STATEMENT  (oi  <#|J«  Rmmmet) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (mi  thm  mOmtrmmi  In  BtmmM  30,  It  di Itmemnt  trmm  RapaN) 


It.  KEY  WORDS  (Cmmttmmm  mm  nvwh  aMt  I#  umiiwr  lOmmHtf  Or  Oimmm  numOme) 

Exponential  auto  regressive  and  moving  average  (EARMA) 
Box  plot 

Variance  Reduction 


to.  ABSTRACT  (CnUrnu  an  nnm  NNa  If  mhiit  an.  IMA  *T  AM*  NiSir) 

T?he  object  of  this  thesis  is  to  examine,  via  simulation, 
the  properties  of  a class  of  single-server,  first-come  first- 
served  queues  in  which  the  service  times  and  interarrival  times, 
while  still  marginally  exponential,  are  auto-correlated  and 
cross-correlated.  The  correlation  is  introduced  by  letting 
the  service  and  interarrival  sequences  be  EARMA-type  processes, 
where  EARMA  stands  for  exponential  autoregressive,  moving  — ° 


oyp 


I JAN  71 


COITION  OP  • NOV  AC  IS  OOCOLRTI 
t/n  0 101*014*  iAOl  | 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  O 


(20.  ABSTRACT  Continued) 


average  sequence.  An  extension  of  these  ideas  brings  in 
the  cross-correlation.  The  waiting  times  in  the  dependent 
queue  are  compared  to  the  waiting  times  (with  known 
distribution)  in  the  independent  M/M/1  queue  using  data 
analytic  and  formal  statistical  methods.  Variance 
reduction  techniques  are  also  studied;  these  use 
distributionally  known  aspects  of  the  M/M/1  queue  (simulated 
with  the  same  exponential  sequences  as  the  dependent  queue) 
to  control  the  unknown  aspects  of  the  dependent  queue. 


f ACC£S$!0‘  ~ I 

NTIS 

: Ciit:on  p* 

DOC 

Section  □ 

UHANt:: 

□ 

jusunc 



«ffKBUT/;;/:.vAiLAMLirr  tog 

Blit.  AVi'.it-  e.-ui/Qf  SPECIAL 


Approved  for  public  release;  distribution  unlimited. 

A Simulation  Study  of  a 
Class  of  First-Come  First-Served 
Queues  with  EARMA  Correlation  Structure 


Prasert  Boonsong 

Lieutenant  Commander,  Royal  Thai  Navy 
B.S.,  Royal  Thai  Naval  Academy,  1964 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  1978 


Author 


Approved  by: 


Thesis  Advisor 


Second  Reader 


Chairman/*  Department  of  Operations  Research 


)ean  dr  Inf ormation? and  Policy  Sciences 


ABSTRACT 


The  object  of  this  thesis  is  to  examine,  via  simulation, 
the  properties  of  a class  of  single-server,  first-come 
first-served  queues  in  which  the  service  times  and  inter- 
arrival times,  while  still  marginally  exponential,  are 
auto-correlated  and  cross-correlated.  The  correlation  is 
introduced  by  letting  the  service  and  interarrival  sequences 
be  EARMA-type  processes,  where  EARMA  stands  for  exponential 
autoregressive,  moving  average  sequence.  An  extension  of 
these  ideas  brings  in  the  cross-correlation.  The  waiting 
times  in  the  dependent  queue  are  compared  to  the  waiting 
times  (with  known  distribution)  in  the  independent  M/M/1 
queue  using  data  analytic  and  formal  statistical  methods. 
Variance  reduction  techniques  are  also  studied;  these  use 
distributionally  known  aspects  of  the  M/M/1  queue  (simulated 
with  the  same  exponential  sequences  as  the  dependent  queue) 
to  control  the  unknown  aspects  of  the  dependent  queue. 
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I.  INTRODUCTION 

The  object  of  this  thesis  is  to  examine,  via  simulation, 
the  properties  of  a class  of  single-server  first-come  first- 
served  queues  in  which  the  service  times  and  interarrival 
times,  while  still  marginally  exponential,  are  auto-correlated 
and  cross-correlated.  The  correlation  is  obtained  by  making 
the  service  and  interarrival  time  sequences  multivariate  EARMA 
processes.  An  ancillary  aspect  of  the  study  is  the  investi- 
gation of  the  efficiency  of  a number  of  control  variable 
variance  reduction  schemes.  These  schemes  are  made  possible 
by  the  simplicity  of  the  structure  of  the  multivariate  EARMA 
processes,  and  the  fact  that  the  null  case  (no  correlation) 
gives  the  M/M/1  queue,  whose  properties  are  known  analytically . 

An  immediate  problem  which  arises  in  the  simulation  is 
that  it  is  not  known  how  far  out  along  the  sample  path  the 
simulation  must  be  carried  in  order  for  the  steady  state 
properties  of  the  queue  to  hold.  This  is  a general  problem 
which  arises  in  queuing  simulations  and,  in  order  to  cope  with 
it,  the  simulations  are  run  in  500  parallel  replications. 

This  allows  one  to  use  many  methods  of  time  series  analysis 
and  data  analysis  to  obtain  graphical  verification  in  a se- 
quential manner  of  the  approach  to  a steady  state.  The  repli- 
cations also  allow  for  examination  of  the  complete  steady 
state  waiting  time  distribution  in  the  queue,  and  not  merely 
the  mean  waiting  time. 


» 


While  all  of  these  aspects  of  the  queuing  simulation 
are  interwoven#  we  can  separate  them  out  for  purposes  of 
discussion  as  follows. 

(i)  An  elementary  aspect  of  the  queue  which  can  be  examined 
is  the  mean  waiting  time.  The  most  efficient  (smallest 
variance)  estimate  of  this  is  the  sample  path  average 
which  is  again  averaged  over  replications.  This  repli- 
cation of  course  reduces  the  variance  by  a factor  of 

m,  where  m is  the  number  of  replications.  Furthermore 
the  existence  of  replications  allows  one  to  look  at  the 
complete  distribution  of  the  sample  path  average  esti- 
mate. This  distribution  should  be  normally  distributed, 
possibly  even  before  a steady  state  has  been  reached. 

The  object  of  this  part  of  the  study  is  to  see  how 
the  mean  steady  state  waiting  time  in  the  queue  varies 
with  the  traffic  intensity,  t,  and  the  correlation 
parameters,  and  in  particular,  how  this  mean  waiting 
time  differs  from  the  waiting  time  in  the  M/M/1  queue 
(null  case) . 

(ii)  A more  detailed  aspect  of  the  queue  which  has  to  be 
examined  is  the  complete  distribution  of  the  steady 
state  waiting  time.  In  the  M/M/1  queue  this  is 
exponential,  given  that  the  waiting  time  is  greater 
than  zero. 

The  first  problem  here  is  to  determine  that  one 
is  actually  looking  at  the  steady  state  waiting  time. 
This  is  done  by  examining  the  empirical  distribution 
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function  of  the  waiting  times  across  replications 

at  times  n^  and  n2  (n2  > n^) , where  the  index  n 

refers  to  the  customer  number.  Of  course  the  two 

samples  from  W and  W , each  of  size  m,  will  be 
nl  n2 

correlated  and  this  makes  standard  statistical  methods 
for  comparing  samples  invalid.  However  the  correla- 
tion between  these  samples  can  be  calculated  in  order 
to  determine  that  n2  is  sufficiently  large  relative  to 
n^  for  the  correlation  to  be  negligible.  In  actual 
fact  the  comparison  is  done  by  looking  at  box  plots 
of  the  W 's  at  successive  values  n.  < n_  < n,  < ...  . 
This  method  is  graphic  and  allows  the  simulator  to 
see  how  the  steady  state  is  being  approached. 

(iii)  A third  problem  which  is  involved  in  looking  at  the 
complete  distribution  of  the  waiting  time  is  that  a 
sample  of  size  m = 500  is  not  sufficient  to  really 
get  a good  estimate  of  the  steady  state  distribution 
of  the  waiting  time.  Thus  successive  samples  along 
the  parallel  sample  paths  are  combined,  where  it  is 
determined  that  the  samples  are  far  enough  apart  to 
be  approximately  independent. 

Again  the  object  here  is  to  see  how  the  complete 
distribution  of  the  steady  state  waiting  time  differs 
from  the  exponential  distribution  (null  M/M/1  case) 
as  traffic  intensity,  t,  and  correlation  parameters 
vary.  In  particular  heavy  traffic  theory  oays  that 
this  distribution  should  be  approximately  exponential 
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as  the  traffic  intensity  parameter  approaches  1. 

The  difficulty  in  verifying  th's,  and  the  heavy 
traffic  approximation  to  the  mean  waiting  time,  is 
that  the  time  to  reach  the  steady  state  becomes  very 
long  as  t -*■  1.  An  additional  question  that  arises 
then  is  whether  it  is  possible  to  start  the  simulation 
in  such  a way  (approximate  steady  state  conditions) 
so  as  to  speed  up  this  convergence.  It  should  be 
remembered,  that  the  queue  with  correlation  is  not 
regenerative  and  no  initial  conditions  for  stationarity 
are  known  or  are  likely  to  be  found. 

(iv)  Variance  reduction  is  always  helpful  in  a simulation 
of  this  kind  since  real  detail  is  always  masked  by 
sampling  fluctuations  i.e.,  the  variability  in  the 
estimate  of  the  unknown  mean  waiting  time.  Because 
of  the  simple  probabilistic-linear  structure  of  the 
EARMA  processes,  running  on  M/M/1  queue  in  parallel 
with  the  correlated  queue  using  common  exponential 
variates  provides  a very  good  control  variable  for  the 
waiting  time  if  the  correlation,  in  a rough  sense, 
in  the  EARMA  queue,  is  small.  However,  as  this 
correlation  increases,  the  amount  of  control  decreases. 
Thus  several  other  control  variables  are  examined  and 
combined  into  a multiple  control  for  the  correlated 
queue.  In  this  way  it  is  hoped  to  obtain  variance 
reduction  for  any  values  of  the  correlation  parameters 


in  the  EARMA  queue.  The  actual  degree  of  control 
which  can  be  attained  is  examined,  using  the  repli- 
cated simulation,  as  a function  of  the  parameters 
of  the  queue. 
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The  first  order  moving  average  process  and  the  first- 
order  autoregressive  processes  were  combined  by  Jacobs  and 
Lewis  (1977)  to  form  the  EARMA(1, 1)  sequence.  We  first 
describe  these  two  first  order  sequences,  and  then  the  method 
of  combining  them. 

A.  THE  EXPONENTIAL  AUTOREGRESSIVE  PROCESS,  EAR(l) 

The  standard  linear,  first-order  autoregressive  model  for 
a stationary  sequence  of  random  variables  {X^}  is  defined 
by  the  equation 

(II. 1)  X±  = pXj,_1  + ei  ; i = 0,  ±1,  ±2,  ...  , 

where  p is  a constant  which  is  less  than  1 in  absolute  value 
and  {e^}  is  a sequence  of  independent  and  identically  dis- 
tributed random  variables.  Gaver  and  Lewis  (1978)  showed 
that  if  the  {X^}  sequence  were  to  have  an  exponential  marginal 
distribution  with  parameter  X , then  the  parameter  p should  be 
greater  than  or  equal  to  zero  and  less  than  one,  and  e^ 
should  be  zero  with  probability  p or  an  exponential  (X) 
random  variables,  E^  with  probability  1-p.  Thus 

Xj.  ” pX^i  + e ^ } 1 * 0,  il,  i2,  ...  , 

becomes 
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i * 0,  ±1,  ±2,  ...  where  {Ei>  is  an  identical  independent 
distribution  (i.i.d.)  sequence  of  exponential  (X)  random 
variables.  Note  that  for  this  EAR (1)  model  the  distribution 
of  the  depend  on  p,  the  multiplicative  weight  of  X^_^. 
Also  is  not  an  absolutely  continuous  random  variable. 

Thus  standard  results  (Mallows,  1968)  to  the  effect  that  the 
X.^  in  an  autoregressive  process  become  approximately  normal 
as  p -*■  1 do  not  hold.  In  fact  in  the  EAR (1)  process  (II. 2), 
the  X^  are  exponentially  distributed  for  0 £ p < 1. 


B.  THE  EXPONENTIAL  MOVING  AVERAGE  PROCESS,  EMA(l) 

The  first  order  moving  average  exponential  process  EMA(l) 
was  developed  by  Lawrance  and  Lewis  (1977)  in  the  form  of 

0Ei  w.p.  0 

(II. 3)  Xi  * 

3E^  + E,^^  w.p.  (1-0)  , 

i=0,  ±1,  ±2,  ...  , where  3 is  greater  than  or  equal  to  zero 
and  less  than  or  equal  to  one,  and  {E^}  is  attain  a sequence 
of  i.i.d.  exponential  (X)  random  variables.  (This  is  the 
backward  case  — a forward  case  is  also  possible.)  The  X^'s 
have  an  exponential  marginal  distribution  and  are  only  serially 
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dependent  for  lag  one;  this  model  is  highly  tractable  and 
a full  account  of  the  statistically  useful  properties  was 
obtained  by  Lawrance  and  Lewis  (1977) . The  forward  model 
combines  E^  with  E^+^  instead  of  with  E^_^. 

C.  THE  EXPONENTIAL  AUTOREGRESSIVE-MOVING  AVERAGE 
PROCESSES  EARMAd/l) 

The  EAR (1)  model  of  Section  IIA  and  the  EMA(l)  model 
of  Section  IIB  are  combined  in  the  following  form  to  give 
an  EARMAd/l)  process.  This  is  a non-Markovian  model  which 
contains  the  EMA(l)  and  EAR(l)  models  as  a special  case.  The 
defining  equations  for  the  EARMA (1, 1)  process  are 

8Eif  w.p.  8, 

(II. 4)  Xt  = 

6Ei  ^ / w.p.  1—8/ 

(0  < 8 ^ 1;  1 ® 0 , il , ±2 / ...) 


where 


pA^_  2 w.p.  p, 

Ai-1  " 

pA^_2  + Ei-i»  w.p.  1-p, 

(0  < p < 1;  i » 0,  ±1,  ±2/  ...) 

and  {E^}  is,  as  usual,  a sequence  of  i.i.d.  exponential  (X) 
random  variables.  Thus  instead  of  E^  being  combined  with  the 
exponential  random  variable  E^_^,  it  is  combined  with  an 
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exponentially  distributed  random  variable  which  is  an  EAR(l) 
combination  of  E^_^,  Ej_2'  Ej__3'  •••  • The  serial  correla- 
tions for  the  EARMA  (1, 1)  process  are  given  in  Jacobs  and 
Lewis  (1977)  as 

p(j)  * corr(X^X^+j)  ■ c(p,6)p^  1 

(j  = 0,±1,±2,±3,  0 < p < 1?  0 < B < 1) . 

where  c(p,8)  * 6(1-6)  (1-p)  + (1-6)2P. 

Note  that  when  p = 0 the  EABMA (1,1)  process  reduces  to  an 
EMA(l)  process;  if  8 = 0 it  is  the  autoregressive  EAR(l) 
process,  and  when  6 3 1 it  is  the  usual  sequence  of 
exponential  (A)  random  variables.  We  will  say  that  Xi 
is  autoregressive  over  E^,  E^_^,  ...  . 


III.  USE  OF  MIXED  STRUCTURES  EAKMA(1,1) 

IN  MODELLING  QUEUES 

• 

Consider  for  simplicity  a queue  with  a single  input 
stream  and  a single  server  and  a first-come-first-served 
(FIFO)  service  discipline.  Let  S^,  i 3 0,  1,  2,  ...  , 

denote  the  service  time  for  the  ith  arrival,  and  let  X^, 

i * 1 , 2 , ...  , denote  times  between  arrival  of  the  ith  and 
(i-l)th  customers  (interarrival  times).  As  is  usual  we  assume 
that  the  first  customer  (with  service  time  SQ)  arrives  at 
time  zero  and  finds  the  queue  empty. 

If  the  (Si>  and  {X^  sequences  are  i.i.d.  exponential 
random  variables  with  parameters  X and  a respectively,  we 
have  M/M/1  queue. 

Now  let 

E^  be  i.i.d.  exponential  (X);  i = 0,  1,  2,  ...  , 

be  i.i.d.  exponential  (a);  i = 0,  1,  2,  ...  . 

We  want  to  model  queues  with  correlated  (autocorrelated 
and/or  cross-correlated)  service  and  inter-arrival  times, 
the  service  and  inter-arrival  times  both  having  marginally 
exponential  distribution  (MD3WD/1  queue) . Exponential 
marginals  is  a common  assumption  and  with  it  the  MDxMD/1 
queuing  model  includes  the  M/M/1  queue  as  a special  case. 
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The  dependence  in  the  arrival  and  service  processes 
is  created  here  by  defining  and  as  a bivariate  depen- 
dent sequence  of  random  variables  with  exponential  marginal 
distributions.  This  is  done  by  letting  {S..}  be  EABMA(1,1) 
over  E^,  (a/X)e^,  (a/X)e^_^,  ...  , i.e. 


so  - Eo 


BE. 


w.p.  B; 


(a/X) e1; 


B*i  + A^  w.p.  (1-8) ; 


BE. 


w.p.  8; 


pA. 


w.p.  p; 


6E2  + A2,  w.p.  ( 1— 6 ) ; 


pA1  + j£2  w.p.  ( 1— p ) ; 


Si 


BEi# 


w.p.  8; 


PA 


i-1 


w.p.  p; 


Ai 


®Ei  + Ai'  w,p*  ? pAi-l  + X®i  w*p* 


We  also  let 


although  one  could  further  make  the  X^'s  correlated.  With 
the  present  assumptions  the  input  process  is  still  a Poisson 
process,  as  in  the  M/M/1  queue.  The  cross-coupling  between 
the  sequence  {S^}  and  {X^>  is  apparent  from  the  structure. 

Note  that  the  S^'s  are  now  correlated  because  of  the  cross- 
coupling. The  {Si>  is  a sequence  of  dependent  exponential 
random  variables,  but  not  an  EARMA(1,1)  sequence.  It  is  a 
type  of  pseudo-EARMA (1,1)  sequence. 

The  interpretation  of  this  queuing  model  is  that  the 

server  tends  to  speed  up  if  the  queue  gets  long  (in  the 

\ 

past) . Of  course  he  also  slows  down  when  the  queue  gets 
short.  In  the  case  where  p * 0 and  0 = 0 then  the  correlation 
between  {S^}  and  {X^}  will  be  equal  to  one,  i.e.,  p 

* (a/A)X^.  Also  when  0*1,  and  for  any  value  of  p, 
the  queuing  model  will  be  identical  to  the  M/M/1  queue.  Thus 
the  M/M/1  is  included  as  a special  case.  Some  more  complicated 
schemes  for  coupling  the  service  and  inter-arrival  processes 
are  given  in  Lewis  and  Shedler  (1978) . Jacobs  (1978a,  1978b) 
has  studied  heavy  traffic  approximations  for  EARMA-type  queues. 

Other  attempts  to  model  M/M/1  (FIFO)  type  queues  with 
correlated  service  and  interarrival  times  have  been  made. 

None  of  them,  however,  have  given  fixed  (exponential) 
marginal  distributions  for  Sj.  and  or  cross-coupling  between 
these  sequences.  There  are,  of  course,  Markovian  schemes 
based  on  birth-death  representations  which  give  FIFO  queues 
in  which  the  service  times  and  interarrival  times  are  correlated 


RMSyj 

n;  • rv 
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(Cox  and  Smith,  1961) , but  the  scheme  is  quite  different 
from  that  given  here.  In  particular  the  interarrival  and 
service  times  are  not  exponentially  distributed,  and  this 
may  not  conform  to  what  is  observed  in  practice. 


IV.  THE  SIMULATION  MODEL 


A.  INTRODUCTION  AND  NOTATION 

Since  both  queues  we  are  considering,  the  M/M/1  queue 
and  the  MDxMD/1 , are  single-server,  first-in-first  out  queues, 
the  waiting  times  Wr  in  both  are  given  by  the  recursion 
equation 


Wn+1  = max(Wn  + Sn  " Xn+1'0)  ; n “ °'1'2'  * 


These  equations  are  used  to  generate  successive  W^’s  in  the 
simulation.  This  is  the  only  aspect  of  the  queue  which  will 
be  considered.  To  have  looked  at  quantities  like  the  number 
of  customers  in  the  queue  at  departure  times  would  have  com- 
plicated the  programming  and  taken  away  from  the  primary  aim 
of  the  thesis,  which  is  to  explore  the  effect  of  correlation 
on  the  queue  and  to  use  some  relatively  modern  statistical 
methodology  to  do  this. 

We  denote  the  waiting  time  of  the  nth  customer  to 
arrive  in  the  jth  realization  of  the  queue  by  Wn(j).  The 
usual  sample  path  estimate  for  the  mean,  E(w),  of  the 
limiting  distribution  of  waiting  times,  Fw(x) , where 


v*> 


lim  P{W_  < x)  , 
ir""  n ~ 


mm 


is  the  sample  path  average 


Wn(j) 


I I*tU> 


which  can  be  shown  to  converge  to  the  mean  E(W)  as  n -►  ®. 

The  M independent  sample  path  realizations  Wn(j), 
n » 0,1,2,  ...,  j - 1,2,  ...,  M can  be  used  to  obtain  estimates 


(1)  the  distribution  of  the  estimate  WnCj).  This  should 
be  normal  for  large  n.  This  can  be  examined  from 
the  sample  Wn(j),  j = 1,  ...,  M.  In  addition  the 
mean  of  this  sample 


S l 


M n 

1ST  l l 

j=l  1=0 


is  a grand  mean  estimating  E (w)  whose  variance  is 
computable  from  the  sample  W^(j),  j s 1, 


Thus 


Var (W£) 


jjj(sample  variance  of  W^(j)'s) 


5<mTT  J 


w/> 


(2)  The  distribution  of  (not  only  its  mean)  from  the 
sample  Wn(j),  j * 1/  ••••  M*  Since  this  is  a random 
sample,  all  the  classical  methods  are  available, 
e.9.,  histograms,  normal  plots,  empirical  c.d.f's. 

(3)  The  correlations  between  successive  waiting  times, 
such  as  corr(WT  (j),Vf  (j)),  nn  < n0.  This  is 

nl  n2  12 

given  as 


B<  l (Mn 
j = l 1 


V*  > <Wn 

nl  n2 


- W ) ) , 
n2 


divided  by  the  sample  standard  deviations  of  the  two 
samples,  where 


wn  = » l (W  (j)) 

n.  « “ n, 

1 3*1  * 

is  am  average  across  replications. 


B.  PROGRAM  STRUCTURE 


1.  General 


To  simulate  this  particular  EARMA(1,1)  with  cross- 
correlated  service  time  amd  inter-arrival  time,  a FORTRAN 
prograun  has  been  written  to  calculate  the  mean  waiting  time 
Wn(j).  It  also  calculates  the  mean  of  various  statistics 
for  the  M/M/1  queue  using  the  same  exponential  deviates  that 
are  used  to  generate  the  correlated  queue.  These  statistics 
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are  used  to  control  the  estimates  from  the  MDxMD/1  queue  for 
purposes  of  variance  reduction.  The  maximum  number  of 
customers  which  the  program  can  handle  for  each  replication 
is  n = 10,000.  This  computation  is  divided  into  a maximum  of 
10  steps  and  the  samples  of  waiting  times  at  these  points  are 
stored  in  a file.  Thus  for  example,  if  we  use  the  maximum 
size  n * 10,000,  we  have  ten  data  sets  w1000(j)f  W2000  (^' 

...,  W^q  000^'  3 . . . , M.  Thus  the  time  evolution  of 

the  queue  cam  be  studied.  Moreover  there  is  a facility  to 

restart  the  simulation  to  get  wn ,000  ^'  W12,000  ^}  ••• 

This  way  the  evolution  of  the  queuing  process  can  be  studied 
as  far  out  as  needed,  e.g.,  to.  n = 270,000,  or  n = 1,000,000  etc. 

I 

Another  program  will  take  care  of  reading  from  the  output  file 
to  do  data  analysis,  i.e.,  box  plot  of  waiting  times  at 
various  steps,  average  waiting  times  W (j)  and  their  statistics. 

2.  MAINl  Program 

a)  The  variables  are  to  be  read  in  as  follows: 


N = Number  of  arrivals  to  be  generated 

(max  is  10,000) 


RS,  RS 

* 

Arrival,  Service  rate. 

BETA,  RHOS 

3 

0,  P 

KN 

= 

Number 

of  the  run 

KS 

= 

1 

if  first  run,  so 

that  n = 

0 

otherwise  (i.e.. 

restart) 

NREP 

3 

Number 

of  replications 

(max  500) 

NSTEP 

= 

Number 

of  steps  in  this 

run 

(maximum  of  10) 


i_  

d — , 
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NPAST 


I EX,  IES 


IBRTA,  IRHOS 


Number  of  arrivals  of  last  run 
(i.e.,  last  run  stopped  at  n * 120,000) 

Seeds  to  generate  exponential  random 
variable  of  Arrival,  Service  times 

Seeds  to  generate  uniform  random 
variables  to  generate  probabilities 
in  the  EARM(1,1)  process  of  moving 
average  and  auto-regressive. 


b)  The  output  will  be  stored  in  two  separated  files. 

i)  For  continuation  of  the  generation  of  the  waiting 

times  of  each  replication  to  the  next  run  (i.e.,  the  restart), 
the  variables  stored  are  as  follow: 


ARIV 


SERV 


SNEW 


WAIT 


WBAR 


Correlated  queue 


= cumulative  inter-arrival  time  (these  are 
the  same  for  both  queues) 

= cumulative  service  time 

= service  time  of  the  last  customer  in  this 


waiting  time  of  the  last  customer  in  this 


ALAST  = 


cumulative  waiting  time 

Auto-rsgressive  component  of  the  inter- 
arrival time  of  the  last  customer  in  this 


TL  =»  Number  of  customers  who  have  arrived  since 
the  last  zero  waiting  time,  i.e.,  backward 
recurrence  time  in  the  renewal  process  of 
times  of  emptiness  in  the  queue 

TR  = Number  of  zero  waiting  times  up  to  the 
termination  point  of  this  run. 

Uncorrelated  queue 

The  variables  are  defined  as  the  same  as  for 
the  correlated  queue;  in  the  program  the  quantities  are 
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suffixed  by  M,  e.g.  WAITM  is  the  waiting  time  for  the 
uncorrelated  queue  and  WAIT  the  waiting  time  for  the 
correlated  queue . 

ii)  For  analyses  and  displays  of  the  data  in 
each  run  and  each  step  the  data,  including  the  data  with 
control  variables  added,  are  stored  in  another  file  in  the 
sequence  of  waiting  time  in  each  replication  of  the  corre- 
lated queue,  average  waiting  time  in  each  replication  of 
the  correlated  queue,  waiting  time  in  each  replication  of  the 
uncorrelated  queue,  average  waiting  time  in  each  replication 
of  M/M/1  queue,  waiting  time  (controlled)  and  average  waiting 
time  (controlled)  of  the  correlated  queue. 

' 3.  Subroutine  CARMA 

This  subroutine  is  used  to  create  a string  of  service 
times,  Sn,  which  are  cross-correlated  to  the  inter-arrival 
time  Xn,  and  also  to  create  the  moving  average  structure. 

4.  Subroutine  CONVAR,  COPYSM,  COPYMV  and  MPROCV 
These  four  subroutines  are  used  to  compute  the 

statistics  from  the  M/M/1  queue  which  are  used  to  control 
the  data  from  MDxMD/1  queue  for  purposes  of  variance  reduction 

5.  MAIN2  Program  with  Subroutine  COMPAR,  FILL  and  STAT 


These  programs  and  subroutines  are  used  to  display 
and  analyze  the  data  for  each  run.  They  produce  the  figures 
given  in  this  thesis.  (See  Appendix  II.) 


a.  Description  of  Box  Plot 

Subroutine  COMPAR  will  produce  vertical  box 
plots  (McNeil  1977)  parallel  to  each  other  and  write  out 
the  minimum  and  the  maximum  value  of  all  the  data  in  the 
array . 

Each  box  plot  is  obtained  by  first  calculating 
the  lower  and  upper  quartiles  and  the  median  of  the  batch 
of  numbers. 

Where  the  median  of  a batch  of  numbers  is  the 
value  for  which  half  of  the  numbers  in  the  batch  are  larger 
and  half  are  smaller,  the  lower  quartile  is  the  value  that 
divides  the  batch  into  two  parts,  with  1/4  of  the  numbers 
below  this  value,  and  3/4  above  it.  Similarly  the  upper 
quartile  is  the  value  for  which  3/4  of  the  numbers  are 
below  it  and  1/4  above  it. 

The  narrow  rectangular  box  with  ends  corres- 
ponding to  the  lower  and  upper  quartiles  contains  the  median 
point  (*) . The  height  of  the  box,  the  interquartile 
distance,  D , is  then  measured  out  on  each  side  of  the 
quartiles.  Then  the  lowest  and  highest  data  points  that 
fall  between  these  quantities  are  also  marked  by  crosses  (X) . 
Finally,  any  numbers  whose  positions  are  outside  these  crosses 
are  marked  with  circles  (0),  those  more  than  1.5  inter- 
quartile distances  outside  each  quartile  getting  the 
symbol  (@).  The  digit  besides  the  "outlier"  positions  is - 
to  indicate  the  number  of  data  points  that  are  too  close 
to  be  marked  separately.  (See  Figure  1 .) 


HIGHEST 
POINT 
IN  BOX 


0 

O 

0 

X 


MEDIAN 


ALL  POINTS 
MARKED  BY  0 


X 

0 

0 

0 

0 

0 

0 

0 

0 


r 

D/2 


D/2  ^ LOWEST 
POINT 
IN  BOX 


1 


ALL  POINTS 
MARKED  BY  0 


FIGURE  1.  Description  of  Box  Plot 
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V.  ANALYSIS  OF  SIMULATION  RESULTS 

A.  INTRODUCTION 


I 

i | 


■ 

i 

i 


I 


The  object  of  the  simulation  was  to  examine  the  effect 
of  the  cross  correlation,  indexed  by  the  coefficients  p and 
B in  the  EARMA(1,1)  process,  on  the  stationary  waiting  time 
W for  various  values  of  the  traffic  intensity  t.  Both  the 
mean  and  distribution  of  W were  compared  to  the  known  mean 
and  distribution  of  the  waiting  time  in  an  M/M/1  queue 
(same  t and  same  mean  service  time) . 

To  perform  this  comparison  we  initially  simulated  this 
model  for  values  .25,  .50,  .75,  .90  of  the  taffic  intensity 
(t) , for  values  0.0,  .25,  .50,  .75,  & .90  of  the  autoregressive 
parameter  (p) , and  the  single  value  0.50  for  the  moving  aver- 
age coefficient  (3)  for  each  combination  of  t and  p.  These 
values  were  chosen  to  limit  the  amount  of  computing  to  a 
physically  feasible  range,  with  the  idea  that  higher  values 
of  p and  t would  be  examined  if  time  permitted.  The  number 
of  customers,  n,  until  steady  state  is  reached  goes  up  drama- 
tically as  p and  t increase.  The  simulations  are  run  in  500 
parallel  replications,  and  the  10,000  arrivals  for  each  run 
are  divided  into  10  steps  of  1,000  arrivals. 

Successive  groups  of  10,000  arrivals  were  run  for  each 
combination  of  p and  t until  it  was  determined  from  an  analy- 
sis of  the  output  that  the  simulation  had  reached  a steady 
state.  This  determination  was  made  from  box-plots  of  samples 
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W (j),  W (j),  j * 1,  . ..,  M,  for  n,  > n, , from  plots  of 

.1  _ n2  * 1 

W and  W , and  from  formal  statistical  comparisons  of  the 
n n 

two  samples  with  ^ large  enough  compared  to  n^  that  the 
correlation  between  the  two  samples  is  negligible. 

B.  ESTIMATION  OF  TIME  TO  STEADY  STATE 

For  each  run,  the  distribution  of  the  sample  path  average 
estimate  Wn(j)  for  each  step  is  observed  by  the  method  of 
box  plots.  By  comparing  these  plots  along  the  n axis  at 
steps  of  arrivals  of  1,000,  convergence  to  steady-state  and 
normality  can  be  determined.  Figures  la-g  show  the  conver- 
gence to  steady  state  of  the  M/M/1  queue  at  t = .75.  For 
the  EARMA- type  MDxMD/1  queue  with  t = .75,  p = .75,  6 = .50, 
the  box  plots  are  shown  in  Fig.  2a-g.  The  number  of  arrivals 
at  the  estimated  steady  state  for  the  M/M/1  queue  and  for 
the  correlated  queue  for  any  combination  of  p and  t are 
shown  in  table  1. 

In  Figures  1 and  2 note  that  the  variance  of  Wn(j)  is 
decreasing  as  1/n,  so  that  there  is  a continual  shrinkage 
from  left  to  right  in  the  plots.  The  symmetry  and  lack  of 
extreme  values  gives  an  indication  of  normality  of  the  esti- 
mates. Note  that  the  scales  on  successive  figures  change 
through  la  to  lg  etc.  This  is  because  the  plots  are  not 
commensurate.  Range  of  values  is  indicated  by  the  values 
at  top  and  bottom  on  the  left,  e.g.,  2.3460522  and  0.3639930 
in  Figure  la. 
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Correlation 


Traffic  intensity,  t 


.25 

o 

in 

• 

in 

r» 

• 

.90 

o 

o 

• 

40,000 

40,000 

50,000 

50,000 

• 

fO 

cn 

40,000 

50,000 

50,000 

70,000 

.50 

40,000 

50,000 

50,000 

70,000 

in 

r- 

• 

40,000 

50,000 

70,000 

100,000 

.90 

60,000 

70,000 

100,000 

100,000 

Table  1.  Number  of  arrivals  that  the  simulation 
run  indicates  are  needed  approximately 
to  reach  steady  state  (for  8*0.5).  These 
are  read  off  of  plots  such  as  those  given 
for  t * 0.75,  p * 0.75  and  8 * 0.5  in 
Figures  2a-g. 
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Figure  la.  Box  plots  derived  from  500  replications,  j«l,...,500  of  the 
sample  path  estimate  ffn(«J)  of  the  wan  waiting  time.  The 
mean  of  these,Wn,  Is  an  overall  estimate  of  the  mean  waiting 
time.  The  medians  of  the  Wn(j)  are  given  by  the  * in  the 
box  plot. 

Uncorrelated  queue  (M/M/1) 

N « 1 ,000  to  10,000  In  steps  of  1 ,000 
t * 0.75,  1/E(S)  ■ Service  rate  ■ 4.0 
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0 • 57  £ C 761 

Figure  1b.  Box  plots  derived  from  500  replications,  j»l,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these, Wn>  is  an  overall  estimate  of  the  mean  waiting 
time.  The  medians  of  the  Wn(j)  are  given  by  the  * in  the 
box  plot. 

Uncorrelated  queue  (M/M/1 ) 

N - 11, tod  tol^OO  In  steps  of  1,000 
t * 0.75,  1/E(S)  * Service  rate  3 4.0 
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Figure  lc.  Box  plots  derived  from  500  replications,  j*1 500  of  the 

sample  path  estjjnate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  VL«  is  an  overall  estimate  of  the  mean 
waiting  time.  Tne  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plot. 

Uncorrelated  queue  (M/M/1 ) 

N * 21 ,000  to  30,000  In  steps  of  1,000 
t * 0.75,  1/E(S)  * Service  rate  * 4.0 
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Figure  Id.  Box  plots  derived  from  500  replications,  j»l,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  Is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
In  the  box  plot. 

Uncorrelated  queue  (M/M/1 ) 

N * 31,000  to  40,000  In  steps  of  1,000 
t * 0.75,  1/E(S)  * Service  rate  * 4.0 
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Figure  le.  Box  plots  derived  from  500  replications,  j»l 500  of  thi 

sample  path  estate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  Is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Mn(J ) «re  given  by  the  * 
In  the  box  plot. 

Uncorrelated  queue  (M/M/1) 

N * 41 ,060  to  bOTfiOO  In  steps  of  1,000 
t « 0.75,  1/E(S)  - Service  rate  « 4.0 
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Figure  If.  Box  plots  derived  from  500  replications,  j*l,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
In  the  box  plot. 

Uncorrelated  queue  (M/M/1) 

N * Si ,000  to  60,000  In  steps  of  1,000 
t ■ 0.75,  1/E(S)  * Service  rate  * 4.0 
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Figure  lg.  Box  plots  derived  from  500  replications,  jsl,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plot. 

Uncorrelated  queue  (M/M/1) 

N =»  61,000  to  70,000  in  steps  of  1 ,000 
t « 0.75,  1/E(S)  - Service  rate  = 4.0 
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Box  plots  derived  from  500  replications,  j*l,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plots. 

Correlated  queue  (MDxMD/1) 

N * 1,000  to  10,000  in  steps  of  1,000 
t « 0.75,  3 3 0.5,  P 3 0.75.  Service  rate  3 4.0 


Figure  2a 


1.2625055 


0 

c 

0.7005695 

Figure  2b.  Box  plots  derived  from  500  replications,  j*l,...,500  of  the 
sample  path  estjjnate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plots. 

Correlated  queue  (MDxMD/1 ) 

N 3 11,000  to  20,000  in  steps  of  1 ,000 
t * 0.75,  3 = 0.5,  p » 0.75.  Service  rate  * 4.0 
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Figure  2c.  Box  plots  derived  from  500  replications,  j=l,...,500  of  the 
sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wp,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plots. 

Correlated  queue  (MOxMD/l ) 

N a 21,000  to  30,000  in  steps  of  1,000 
t a 0.75,  0 a 0.5,  p a 0.75.  Service  rate  a 4.0 
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Figure  2d.  Box  plots  derived  from  500  replications,  j*l 500  of  the 

sample  path  estimate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
In  the  box  plots. 

Correlated  queue  (MM1D/1 ) 

N = 31,0(50  to  40,000  In  steps  of  1,000 
t ■ 0.75,  3 * 0.5,  p ■ 0.75.  Service  rate  * 4.0 
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0. 186i617 

Figure  2e.  Box  plots  derived  from  500  replications,  j»1 500  of  the 

sample  path  estjjnate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plots. 

Correlated  queue  (MDxMD/1) 

N = 4l .Odd  to  50,000  in  steps  of  1,000 
t * 0.75,  6 - 0.5,  p * 0.75.  Service  rate  • 4.0 
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a 


0 .7982421 

Figure  2f.  Box  plots  derived  from  500  replications,  j*l,...,500  of  the 
sample  path  estjjnate  Wn(j)  of  the  mean  waiting  time-  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
In  the  box  plots. 

Correlated  queue  (MDXMD/I ) 

N ■ 51 ,000  to  60,000  In  steps  of  1,000 
t * 0.75,  8 - 0.5,  p ■ 0.75.  Service  rate  « 4.0  , 
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Figure  2g.  Box  plots  derived  from  500  replications,  j*l 500  of  the 

sample  path  est^nate  Wn(j)  of  the  mean  waiting  time.  The 
mean  of  these,  Wn,  is  an  overall  estimate  of  the  mean 
waiting  time.  The  medians  of  the  Wn(j)  are  given  by  the  * 
in  the  box  plots. 

Correlated  queue  (MD*1D/1) 

N » 61,000  to  70,000  *n  steps  of  1,000 
t ■ 0.75,  8 ■ 0.5,  p ■ 0.75.  Service  rate  * 4.0 
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C.  ESTIMATED  MEAN  WAITING  TIMES 

To  compare  the  change  of  the  estimate  steady  state  mean 
waiting  time  from  the  M/M/1  queue  as  p is  changed,  Table  2a 
shows  the  value  of  mean  waiting  time,  Wn  and  the  estimate 
standard  deviation  of  those  means.  The  mean  service  rate 
is  fixed  at  4.0. 

The  mean  waiting  times  of  the  correlated  queue  for  any 
fixed  value  of  p and  t,  are  plotted  in  Figures  3a  and  3b. 
respectively.  These  show  that  the  mean  waiting  time  increases 

dramatically  as  t + 1,  just  as  for  the  M/M/1  queue  ({5  ■ 1.0), 

* 

for  which  the  mean  waiting  time  increases  as  l/(l-t).  A 
similar  drastic  increase  in  the  mean  waiting  time  because  of 
the  increase  in  cross-correlation,  measured  by  p is  apparent 
in  Figure  3b,  when  t is  large.  Note  that  when  p is  close 
to  one  there  is  very  long  term  correlation  in  the  queuing 
process;  however,  p = 1 is  not  allowed  because  the  system  is 
not  ergodic. 

The  ratio  of  the  estimated  mean  waiting  time,  Wr  for 
the  correlated  queue  to  that  for  the  M/M/1  queue  which  are 


tabulated  in  Table  2b,  are  also  plotted  in  Figure  4a  and 
Figure  4b.  In  particular.  Figure  4a  shows  that  the  correla- 


r | 


p/(l-P)  and  t/(i-t)  so  as  to  be  able  to  predict  E(W)  for 
large  values  of  t.  The  form  of  the  independent  variables 
were. suggested  by  the  fact  that  a log  transform  showed 
lnE(W)  to  go  up  approximately  as  ln(l-p)  for  large  p,  and 
the  fact  that  for  small  p,  E(W)  is  smaller  than  E(W)  in 
the  M/M/1  case.  Also  the  form  t/ (1-t)  was  suggested  by  the 
M/M/1  theory. 

Only  one  value  of  8 was  used  in  the  simulations,  B 3 0.5, 
so  that  the  results  did  not  reflect  this  variable.  The  rough 

fit  obtained  for  t close  to  one  was,  for  the  MDxMD/1  queue 
with  parameters  t,p,E(S). 

(V . 1)  E*  (W)  = j—  t0. 5 + £j§|£]E(S)  . 


One  use  for  this  formula  would  be  to  predict  E* (W)  for 
high  t so  that  a simulation  could  be  started  with  WQ  exponen- 
tial with  mean  E(W  ) (the  probability  of  a zero  waiting  time 

n 

can  be  neglected  for  high  t) . This  reduces  the  transient 
effect,  as  we  will  see  later. 

Formula  (V.I)  predicts  that  for  t = 0.99,  p 3 0.25 
and  E (S)  3 0.25,  the  mean  waiting  time  would  be 


E*  (W) 


14.43  , t 3 0.99,  p * 0.25,  E(S)  3 0.25 


The  mean  for  the  M/M/1  queue  is  given  as 
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E (W) 


24.75 


A simulation  was  performed  for  this  traffic  intensity 
(t  **  .99),  starting  with  WQ  = 0;  the  simulation  was  found  to 
converge  to  steady-state  for  n * 280,000  and  then  (see 
Figures  12a  and  12b) 


280,000 


14.194 


with 


s.d. (W280/ qoq)  0.1731. 


This  is  in  reasonable  agreement. 

The'  effect  of  using  this  value  for  E(WQ)  in  the  simulation 
is  discussed  in  Section  V.D. , 


Correlation 


Traffic  Intensity,  t 


.25 

o 

in 

• 

.75 

.90 

o 

o 

• 

0.073763 

(0.000064) 

0.192998 

(0.000156) 

0.484147 

(0.000501) 

1.258926 

(0.003004) 

.25 

0.077745 

(0.000730) 

0.210338 

(0.000165) 

0.543177 

(0.000597) 

1.444641 

(0.003373) 

.50 

0.083026 

(0.000091) 

0.237792 

(0.000230) 

0.649296 

(0.000899) 

1.798944 

(0.004716) 

.75 

0.091539 

(0.000130) 

0.296498 

(0.000440) 

0.926665 

(0.001709) 

2.819429 

(0.009518) 

.90 

0.101673 

(0.000188) 

0.413239 

(0.000950) 

1.654127 

(0.004754) 

5.717662 

(0.032587) 

M/M/l 

Queue 

0.083365 

(0.000076) 

0.250211 

(0.000231) 

0.749701 

(0.000983) 

2.254441 

(0.007270) 

Table  2a.  Estimated  mean  W and  standard 
deviation  of  thenestimated  mean 
of  the  average,  steady  state 
waiting  time  in  the  correlated 
queue  for  various  p and  t (6  » 0.5 
and_Ue  * 0.25).  The  estimated  s.d. 
of  wnSis  given  in  brackets. 
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(x  Mean  Service  Time,  E(S) ) 


Traffic  Intensity,  t 

m 

.25 

.50 

.75 

.90 

o 

o 

• 

0.8852 

0.7720 

0.6455 

0.5595 

.25 

0.9329 

0.8414 

0.7242 

0.6421 

.50 

0.9963 

0.9512 

0.8657 

0.7995 

.75 

1.0985 

1.1860 

1.2356 

1.2531 

.90 

1.2201 

1.6530 

2.2055 

2.5412 

Table  2b.  Ratio  of  the  Estimated  mean  Wn  of 
the  average  waiting  time  in  the 
correlated  queue  to  the  known 
average  waiting  time  of  the 
uncorrelated  (M/M/1)  queue  (B  = 0 
for  various  values  of  p and  t. 


Traffic  Intensity  t 


Figure  4a.  Smoothed,  estimated  mean  waiting  time  in  the  correlated 
MDxMD/1  queue  divided  by  the  known  mean  waiting  time  in 
the  M/M/1  queue  as^  a function  of  traffic  intensity  t. 
Thus  the  result  Ts  Tn  uni ts  of- the  mean  waiting  time  in 
M/M/1  queue.  Here  the  correlation  parameter  p is  the 
parameter.  Fixed  8 * 0.5 


2.751 


^ 1.25 


[ill] 


correlation  p 

Figure  4b.  Smoothed,  estimated  mean  waiting  time  in  the  correlated 
MDxMD/1  queue  divided  by  the  known  mean  waiting  time  in 
the  M/M/1  queue  as  a function  of  the  correlation  parameter  p 
The  result  is  in  units  of  the  mean  waiting  time  in  M/M/1 
queue.  Here  the  traffic  Intensity  t Is  the  parameter. 

Fixed  B * 0.5. 


D.  ESTIMATED  DISTRIBUTION  OF  WAITING  TIME  W 

The  simulation  program  was  set  up  to  enable  us  to 

\ 

look  at  statistics  on  W and  W across  the  replications 

n n 

at  various  n's.  In  addition  the  file  system  enabled  us 
to  look  at  joing  properties  of  the  samples  at  two  different 
n's,  and  in  particular  at  correlations.  The  reason  for 
this  is  that  a sample  of  size  500  is  rather  small  when  one 
wants  to  look  at  the  steady  state  distribution  of  Wn, 
namely  W.  One  could  increase  the  number  of  replications, 
but  this  is  constrained  by  the  size  of  the  computer. 

However,  since  the  simulation  has  to  be  carried  out  well 
beyond  the  set  in  of  the  steady  state  for  safety  sake, 
one  might  as  well  pool  samples  of  waiting  times  along 
the  sample  path  if  they  are  sufficiently  far  apart  to  be 
uncorrelated,  or  approximately  so.  Then  instead  of  looking 
at  a sample  size  of  50Q,  we  may  pool  waiting  times  at  10 
steps  of  n,  n1#  n2,  ...,n1Q,  to  make  a sample  size  of  5000. 
Then  the  waiting  time  distribution  for  a particular  run  can 
be  shown  as  in  fig.  5a,  5b;  these  are  actually  histograms 
showing  many  of  the  estimated  sample  statistics  which 
might  be  of  interest. 

• In  each  step  of  the  simulation  runs  which  were  performed 
the  correlation  between  waiting  times  1000  arrivals  apart 
were  obtained  and  found  to  be  negligible,  (see  Table  3a,  3b) , 
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n =*  no.  of 
arrivals 


1.0  -0.0118  0.0 

-0.0118  1.0  -0.0' 
0.0135  -0.0450  1.0 


10000 


0.0135 

0.0002 

-0.0384 

-0.0450 

0.0403 

-0.1079 

1.0 

0.0012 

-0.0088 

0.0012 

1.0 

-0.0491 

-0.0088 

-0.0491 

1.0 

Table  3a.  Correlation  Matrix  of  Waiting  Times  for 

Different  Indices  n,  in  the  Correlated  Queue, 
t = 0.75,  8 « 0.50,  p = 0.75. 


warn* 


n * no.  of 
arrivals 

6000 

7000 

8000 

9000 

10000 

6000 

1.0 

-0.0371 

0.0280 

0.0105 

-0.0131 

7000 

-0.0371 

1.0 

-0.0413 

0.0779 

-0.1707 

8000 

0.0280 

-0.0413 

1.0 

0.0189 

-0.0241 

9000 

0.0105 

0.0779 

0.0189 

1.0 

-0.0977 

10000 

-0.0131 

-0.1707 

-0.0241 

-0.0977 

1.0 

Table  3b.  Correlation  Matrix  of  Positive  Waiting  Times 


(No  Zeros)  for  Different  Indices  n in  the 
Correlated  Queue,  t ■ .75,  8 = -50,  p * 0.75. 


. 


1 


thus  justifying  the  pooled  samples.  Of  course  a more  effi- 
cient use  of  the  data  would  have  been  to  calculate  a correlo- 
gram  of  estimated  serial  correlations  along  each  (stationary) 
sample  path  and  average  over  the  five  hundred  replications 
to  obtain  a more  precise  estimate  of  each  serial  correlation 
together  with  a variance  estimate  of  that  estimate. 

Another  reason  for  looking  at  the  samples  Wn(j ),  j 3 1,  ..., 

500  at  different  values  of  n is  to  assess  the  convergence  to 

a steady  state  in  terms  of  the  whole  distribution  of  Wn, 

rather  than  just  its  mean  value.  There  are  a number  of  ways 

of  doing  this  display  and  a particularly  simple  and  graphic 

way  is  to  give  box-plots  for  the  successive  samples  on  the 

same  scale.  Unfortunately , this  proved  to  be  difficult  with 

presently  available  software;;  thus  in  Figs.  6a-g  we  give 

box  plots  of  Wn(j)  for  n 3 1,000(1,000)70,000.  The  scales 

in  these  figures  are  unfortunately  not  commensurate,  but  the 

convergence  to  steady  state  can  be  seen.  Note  that  unlike 

the  W (j)'s,  the  variance  of  the  W (j)'s  are  not  decreasing 
n n 

with  n;  in  fact  that  should  converge  to  var (W) . If  some  of 
the  distributional  detail  in  the  plots  is  overwhelming,  one 
can  look  solely  at  the  *'s  in 'the  boxes,  which  are  the  median 
values  in  the  successive  samples.  Since  these  samples  are, 
as  we  have  seen,  approximately  uncorrelated,  the  *'s  could 
have  been  smoothed. 

The  plots  in  Figures  6a-g  are  tied  to  zero  by  the  occurrence 
of  zero  waiting  times  (on  the  average  (l-t)%  of  the  sample) . 

Thus  plots  of  the  positive  waiting  times  are  given  in  Figures  7a-g. 
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Figure  5b. 


Histogram  and  sample  statistics  for  combined  waiting  time 
data  (with  2,293  zero  waiting  times  removed)  In  the 
simulated,  uncorrelated  M/M/1  queue,  l.e.  all  Wn(j)'s  for 
j-1 . . ,500  and  n » 51,000(1000)60,000.  t » 0.75,  1/E(S)  * 4.0 
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WAITING  TIME  WITHOUT  ZERC 
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Figure  6a.  Box  plots  of  the  sample  waiting  times  Wn(j),  j*l,...,500 
for  values  of  n * 1,000  to  10,000  in  steps  of  1,000. 
Correlated  queue;  t * 0.75,  6*  0.5,  p*  0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and  minimum 
sample  values. 
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Figure  6c.  Box  plots  of  the  sample  waiting  times  Wn(j),  j«l,...,500 
for  values  of  n * 21,000  to  30,000  In  steps  of  1,000. 

Correlated  queue;  t * 0.75,  6 = 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and  minimum 
/ sample  values. 
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Figure  6d.  Box  plots  of  sample  waiting  times  Wn(j),  j«l 500 

for  values  of  n * 31,000  to  40,000  in  steps  of  1,000. 
Correlated  queue;  t a 0.75,  B * 0.5,  p * 0.75. 
figures  at  top  and  bottom  at  left  are  maximum  and  minimum 
sample  values. 
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Figure  6e.  Box  plots  of  sample  waiting  times  Wn(j),  j*l 500 

for  values  of  n « 41,000  to  50,000  in  steps  of  1,000. 
Correlated  queue;  t a 0.75,  8 * 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and  minimum 
sample  values. 
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Figure  6f.  Box  plots  of  sample  waiting  times  WQ( j ) , j»l 500 

for  values  of  n * 51,000  to  60,000  in  steps  of  1,000. 
Correlated  queue;  t a 0.75,  8 * 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 
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Figure  6g.  Box  plots  of  sample  waiting  times  W (j),  j*l,...,500 
for  values  of  n = 61,000  to  70,000  in  steps  of  1,000. 
Correlated  queue;  t 3 0.75,  0 3 0.5,  p * 0.75. 
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Figure  7a.  Box  plots  of  the  sample  waiting  times  Wn(j),  j»l,...,500 
with  zero  waiting  times  removed  for  values  of 
n * 1,000  to  10,000  in  steps  of  1,000. 

Correlated  queue;  t * 0.75,  6 3 0.5,  p 3 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 


Box  plots  of  the  sample  waiting  times  Wn(j),  j=l,...,500 
with  zero  waiting  times  removed  for  values  of 
n * 1 1 ,d)00  to  20,000  in  steps  of  1 ,000. 

Correlated  queue;  t » 0.75,  B = 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 
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Figure  7c.  Box  plots  of  the  sample  waiting  times  Wn(j),  j*l 
with  zero  waiting  times  removed  for  values  of 

steps  of  1,000. 

Correlated  queue;  t * 0.75,  6 * 0.5,  p ■ 0.75. 
Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 
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Figure  7d.  Box  plots  of  the  sample  waiting  times  Wn(j),  j»l,...,500 
wi th  zero  waiting  times  removed  for  values  of 
n « 3T7000  to  40,000  in  steps  of  1,000. 

Correlated  queue;  t * 0.75,  6 * 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 
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Figure  7e.  Box  plots  of  the  sample  waiting  times  Wn(j),  j*l,...,500 
with  zero  waiting  times  removed  for  val ues  of 
n » 4T7W0  to  56,000  In  steps  of  1,000. 

Correlated  queue;  t - 0.75,  B * 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum  and 
minimum  sample  values. 
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Figure  7f. 


Box  plots  of  the  sample  waiting  times  Wn(j),  j«l,. 
with  zero  waiting  times  removed  for  values  of 
n - 51,000  to  60,000  in  steps  of  1,000. 

Correlated  queue;  t * 0.75,  6 » 0.5,  p » 0.75. 
Figures  at  top  and  bottom  at  left  are  maximum 
and  minimum  sample  values. 
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Figure  7g.  Box  plots  of  the  sample  waiting  times  Wn(j),  j=l,...,500 
wl th  zero  waiting  times  removed  for  values  of 
n » 61,000  to  70,000  in  steps  of  1,000. 

Correlated  queue;  t * 0.75,  B 3 0.5,  p * 0.75. 

Figures  at  top  and  bottom  at  left  are  maximum 
and  minimum  sample  values. 
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Another  question  of  interest  is  how  the  parameters  p and 
t affect  the  waiting  time  in  the  MDxMD/1  queue,  and  how 
these  distributions  differ  from  the  exponential  distribution 
of  the  positive  values  of  for  the  uncorrelated  M/M/1 
queue.  There  are  several  graphical  techniques  for  doing  this, 
and  some  of  the  results  are  as  follows. 

a)  Box  plots  are  used  to  compare  the  distribution  of 
waiting  times  at  8 * 0.5  for  several  fixed  values  of  traffic 
intensity  t as  the  correlation  parameter  p is  varied.  These 
are  shown  in  fig.  8a-d  and  complement  Figure  3b  also  for 
fixed  values  of  the  correlation  parameter  p,  the  variation 
of  the  distribution  of  W is  shown  in  Figs.  9a-f.  Note  that 
in  these  figures  10  approximately  uncorrelated  samples  in  the 
steady  state  have  been  combined. 

b)  The  box  plots  are  not  useful  for  examining  departures 

/ 

from  exponentiality , and  for  this  purpose  the  histograms  and 
sample  statistics  from  HIST  F were  utilized.  For  the  grid 
of  values  of  p and  t discussed  in  (V.B)  the  coefficients  of 
variation,  skewness  and  kurtosis  were  extracted  from  the 
HISTF  outputs  (no  zeros) , and  these  are  compared  to  the 
exponential  values  C(k)  = 1.0,  y^  * 2.00,  y2  * 6.00  in  table  4. 

It  is  clear  that  for  small  p & t the  distribution  is  less 
skewed  than  the  exponential,  but  as  p and  t increase  the 
distribution  becomes  positively  skewed  and  is  certainly  non- 
exponential. Thus  for  p * 0.9  and  t » 0.9  we  have  C(x)  * 1.1325, 
Y^  * 2.3289,  y2  = 3.6449.  One  explanation  for  this  is  that  for 
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Traffic  Intensit 


§ 0.50 


P1  - -0.0476 
p2  » -0.2442 

P3  = -2.0026 


px  * -0.0433 

p2  = -0.2028 

p3  - -1.7623 


px  - -0.0132 
p2  = -0.1963 
p3  « -1.8133 


.1625 

.8135 

.6253 


0.50 

0.75 

-0.1552 

-0.2045 

-0.5133 

-0.5062 

-2.8502 

-1.9100 

-0.1219 

-0.1743 

-0.4925 

-0.5689 

-3.0085 

-3.0471 

-0.0303 

-0.1117 

-0.1725 

-0.3894 

-1.5305 

-2.0960 

+0.1191 

+0.0605 

+0.3999 

-0.0189 

+1.9612 

-0.8715 

40.5023 

40.3245 

+1.6374 

40.7767 

+12.8466 

+5.8286 

-0.1273 

+0.2455 

+4.6589 


0.1448 

0.3550 

2.1357 


-0.0839 

-0.1937 

-1.3652 


+0.0364 

+0.2571 

+2.3012 


40.1325 

+0.3289 

+2.6449 


Pj^  = coefficient  of  variation  - 1.0 


-0.0603 

-0.3000 

-2.223 


P2  = skewness  - 2.0 
p„  = kurtosis  - 6.0 


Table  4.  Deviations  of  the  Estimated  Sample  Coefficients 
of  Variation,  Skewness  and  Kurtosis  for  the 
Correlated  Queue  Waiting  Times  (Without  Zeros) 
from  the  Known  (Exponential)  Values  for  the 
M/M/1  Queue.  For  various  p and  t values. 


76 


p small,  the  service  time  is  decreased  by  the  correlation 
when  many  people  arrive.  The  distribution  is  therefore 
tighter  than  for  the  M/M/1  case.  However  as  p increases  the 
service  times  become  highly  correlated  through  the  cross- 
coupling, and  this  effect,  which  increases  the  mean  waiting 
time  and  the  skewness  of  the  distribution,  is  dominant. 

c)  Another  question  of  interest  is  whether  for  large  t 
the  heavy  traffic  approximation  holds  (Jacobs,  1978)  and 
the  distribution  of  W is  approximately  exponential.  This  is 
difficult  to  examine  because  the  time  for  the  queue  to  reach 
steady  state  is  inordinately  high  for  large  t.  Thus  one  case, 
t = 0.995,  p = 0.25  was  studied  and  the  HISTF  outputs  are 
shown  in  Figures  10a  & 10b.  The  plots  and  the  sample  statis- 
tics still  show  a large  amount  of  underdispersion  relative 

to  an  exponential  distribution.  Thus  convergence  to  an 
exponential  distribution  as  t + 1,  if  it  occurs,  must  be  very 
slow. 

Plots  for  an  extreme  case,  p = 0,  0*0,  are  shown 
for  various  values  of  t from  0.25  to  0.99  in  Figures  lla-e. 

d)  The  question  of  judging  when  the  transient  in  the 
simulation  has  died  out  is  very  complex  and  difficult,  just  like  any 
question  of  detecting  a signal  in  noise.  In  particular  it  is 

to  be  stressed  that  the  more  graphic  output  one  has,  the  better 
off  one  is.  This  is  another  reason  for  looking  at  the  WR(j) 
samples  as  well  as  the  Wn(j)  samples. 

A 

In  Figure  12a  W is  shown  for  the  MD  x MD/1  queue 
n 

with  t » 0.99,  p * 0.25.  The  judgement  was  made,  and  published 


in  Table  1/  that  equilibrium  was  reached  by  the  time  of  the 
40,000th  customer  and  this  appears  reasonable  from  the 
figure,  which  goes  out  to  n » 280,000.  The  plot  of  Wn 
is  given  in  Fig.  12b?  superposing  it  on  Fig.  12a  is  a help  to 
visualization.  It  of  course  has  a larger  transient  than  Wr 
because  it  drags  in  all  the  previous  annual  times,  so  that  the 

A 

plot  of  Wn  gives  a better  picture  of  the  real  steady  state. 

Of  course,  as  in  any  real  simulation,  E(W)  is  unknown  and 

the  asymptote  in  Fig.  12b  is  a help  in  showing  the  convergence. 

Note  too  that  the  normality  in  Wn  decreases  with  n,  but  that 

of  Vf  does  not. 
n 

Note  too  that  since  we  have  determined  that  W 's 

n 

1,000  arrivals  apart  are  almost  uncorrelated,  local  smoothing 

/s 

could  be  applied  to  the  Wn  plot  to  obtain  a better  picture. 

Output  from  a similar  simulation  for  an  M/M/1  queue, 
again  with  t = 0.99  and  500  replications,  is  shown  in  Figs.  12c 
& d.  The  mean  of  W is  known  to  be  E(W)  * 24.75.  However  the 

A 

plot  of  Wn  is  well  above  this  from  approximately  n = 60,000 

A 

to  n = 260,000,  and  since  the  correlation  between  Vf  values 
is  not  very  extensive,  this  probably  indicates  that  the 
transient  is  such  that  E(V»n)  does  not  go  up  monotonically  to 
E (W) , as  one  might  expect  intuitively  or  from  the  simulation 
out  to  n = 240,000,  but  actually  rises  above  E(W)  and  then 
declines.  It  may  even  have  a damped  oscillation.  The  point 
is  that  all  the  data  should  be  examined-  formal  tests  of  the 
similarity  of  the  WR  samples  at  say  W10Q  Q00  and  wno,000  would 
probably  have  indicated  convergence  (i.e.,  similar  distributions). 


Another  overlaying  the  curve  of  on  that  of  Wn 

is  a help  usually,  as  is  the  horizontal  asymptotic  in  Fig.  '12d. 

Note  that  the  initial  transient  in  the  M/M/1  is  longer  tham 

it  is  in  the  MDxMD/1  queue. 

Another  question  which  arises  is  whether  one  can 

start  the  simulation  in  such  a fashion  as  to  reduce  the  length 

of  the  transient.  An  attempt  to  do  this  is  shown  in  Figs.  13a-b. 

Here  am  extreme  case,  t =»  0.995,  was  taken  amd  the  initial 

value  W was  not  taken  to  be  zero,  but  an  exponential  random 
o 

variable  with  mean  29.0.  The  value  29.0  comes  from  the  pre- 
diction formula  IV. 1.  Roughly  speaking  there  appears  to  be 
faster  convergence.  The  simulation  starting  from  WQ  = 0 
takes  a long  time  to  converge  and  was  too  extensive  to 
be  hamdled  on  the  computer. 
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Figure  8a.  Correlated  queue.  Box  plots  of  a sample  of  waiting  times 

JOTTTFl 500;  i*l ,. . . ,10;  with  zero  waiting  times 

removed,  for  t * 0.25  and  five  values  of  p.  For  each  (t,p) 
pair  ni  is  chosen  to  be  large  enough  for  the  queue  to  be  in 
steady  state,  and  10  approximately  uncorrelated  samples  at 


n-j  9 • • • § 


n-j0  are  combined. 


Here  n1  * 31,000,  31,000,  31,000,  31,000,  51,000 

as  p *.0.0,  0.25,  0.50,  0.75,  0.90,  respectively.  8 * 0.5. 
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Figure  8b.  Correlated  queue.  Box  plots  of  a sample  of  waiting  times 

Wn-j(j).  J*1 500;  1=*1 , . . . ,10;  with  zero  waiting  times 

removed,  for  t 3 0.50  and  five  values  of  p.  For  each  (t,p) 
pair  n-j  is  chosen  to  be  large  enough  for  the  queue  to  be  in 
steady  state,  and  10  approximately  uncorrelated  samples  at 
n^ng n1Q  are  combined. 


1-,  * 31,000,  41,000,  41,000,  41  ,000,  61  ,000 
s 0.0,  0.25,  0.50,  0.75,  0.90,  respectively. 
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Figure  8c. 


Correlated  queue.  Box  plots  of  a sample  of  waiting  times 
WnjU)»  3*1,..., 500;  1*1,. ..,10;  with  zero  waiting  times 
removed,  for  t * 0.75  and  five  values  of  p.  For  each  (t,p) 
pair  ni  Is  chosen  to  be  large  enough  for  the  queue  to  be  in 
steady  state,  and  10  approximately  uncorrelated  samples  at 


n^.n2....,n1A  are  combined. 
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j • 0.0.  0.25,  0.50,  0.75,  0.90,  respectively.  6 * 0.5. 
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Correlated  queue.  Box  plots  of  a sample  of  waiting  times 
Wn.(j),  j3! ,. . • ,500;  1-1 , ... ,10;  with  zero  waiting  times 
removed,  for  p * 0.0  and  four  values  of  t.  For  each  (t,p) 
pair  n-j  is  chosen  to  be  large  enough  for  the  queue  to  be 
in  steady  state,  and  10  approximately  uncorrelated  samples 
at  n-j.ng n1Q  are  combined.  0 * 0.5. 
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Correlated  queue.  Box  plots  of  a sample  of  waiting  times 

.BnTOTTFl 500;  1*1,. ..,10;  with  zero  waiting  times 

removed,  for  p * 0.25  and  four  values  of  t.  For  each  (t,p) 
pair  n-|  Is  chosen  to  be  large  enough  for  the  queue  to  be 
in  steady  state,  and  10  approximately  uncorrelated  samples 

_ 4.  _ — _ _ L i I a A P 


at  n-j  .ng,.  • • » n ^ g 


are  combined.  3 * 0.5. 
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Correlated  queue.  Box  plots  of  a sample  of  waiting  times 

WniU),  j=l ,. . . ,500;  i*l 10;  with  zero  waiting  times 

removed,  for  p * 0.50  and  four  values  of  t.  For  each  (t,p) 
pair  n-i  is  chosen  to  be  large  enough  for  the  queue  to  be 
in  steady  state,  and  10  approximately  uncorrelated  samples 
at  n-|  .ng,. . . ,n,n  are  combined.  $ * 0.5. 
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Figure  9e.  Correlated  queue.  Box  plots  of  a sample  of  waiting  times 

Wn.(j),  j=l 500;  i*l 10;  with  zero  waiting  times 

removed,  for  p » 0.90  and  four  values  of  t.  For  each  (t,p) 
pair  ni  is  chosen  to  be  large  enough  for  the  queue  to  be 
in  steady  state,  and  10  approximately  uncorrelated  samples 
at  n^.ng n,ft  are  combined.  $ * 0.5 
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Figure  9f.  Correlated  queue,  p 3 0.0,  8 3 0.0  (service  time  equals 
t times  previous  Inter-arrival  time).  Box  plots  of  a 

sample  of  waiting  times  Wni(j),  j*l 500;  i»l 10; 

with  zero  waiting  times  removed.  Five  values  of  t.  For 
each  t,  n-j  is  chosen  to  be  large  enough  for  the  queue  to 
be  in  steady  state,  and  10  approximately  uncorrelated 
samples  at  n, ,n?,. . . ,n,Q  are  combined. 
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Oa.  Histogram  and  sample  statistics  for  combined  waiting-time 
data  in  the  simulated,  correlated  MDXMD/1  queue; 
t » 0.995,  6 * 0.5,  p a 0.25,  1/E(S)  = 4.0, 
i.e.  all  Wn(j)‘s  for  j=l 500  and  n * 91,000(1000)100,000. 
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Figure  10b.  Histogram  and  sample  statistics  for  combined  waiting-time 
data  ith  18  zero  waiting  times  removed)  in  the  simulated, 
correlated  MDXMD/1  queue;  t * 0.995,  6 a 0.5,  p = 0.25, 
1/E(S)  * 4.0,  i.e.  all  W„(j)'s  for  j * 1.....500  and 
n = 91,000(1000)100,000. 
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Figure  11a.  Histogram  and  sample  statistics  for  combined  waiting-time 
data  (with  zero  waiting  times  removed)  In  the  simulated, 
correlated  MDxMD/1  queue;  t * 0.25,  6 * 0.0,  p 58  0.0, 

1/E(S)  - 4.0,  i.e.  all  WB(j)'s  for  j = 1 500  and 

n » 11,000(1000)20,000.  n 
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Figure  11b.  Histogram  and  sample  statistics  for  combined  waiting-time 
data  (with  zero  waiting  times  removed)  in  the  simulated, 
correlated  MDxMD/1  queue;  t = 0.50,  8 a 0.0,  p * 0.0, 

1/E(S)  * 4.0,  i.e.  all  Wn(j)‘s  for  j ■ 1 500  and 

n * 21,000(1000)30,000. 
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Figure  11c. 


Histogram  and  sample  statistics  for  combined  waiting-time 
data  (with  zero  waiting  times  removed)  in  the  simulated, 
correlated  MDxMD/1  queue;  t * 0.75,  0 = 0.0,  p * 0.0, 
1/E(S)  = 4.0,  i.e.  all  Wn(j)'s  for  j = 1,...,500  and 
n = 21,000(1000)30,000.  n 


9*4 


FFEdIKC  ICS  SAMPLE  SIZE  - 45C2 

204  SSI  AS]  tit  tCC  SSA  45S  435  301  214  141  99  76  47  37  22  1 


i 


LjiAa.xxa.icx4.XAX«-x«.«.&Aa.xa.xxi.JLi4.x«.ia.4.4.«.A 


uouuoou 

I I I I I I 


m*A»<iP4«r»^v 


4 C UJUHUUJUJ 
U -J-O-XO-o 


at  4(i4i(iu^x 

H ■«  "4 

v»  «umuuiuA 


• 

VI 

ivtvuuiyiv 

r-A 

p- 

oouyvjo 

z 

ii  ii 

U 1 

uiujujumjAAi 

X 

V_A 

tft 

X 

UU'U\U«*"0 

*»• 

— • 

-O 

to  tr  co  r—'CW 

o 

I-A 

«4imuvAMUi 

• 

at 

• ••••• 

o 

p- 

rO-O-yrvKM-o 

07 

Z 

O 

« 

Ui 

at 

o 

UI 

• 

«*> 

•**1 

U. 

*•  • 

a 

mv» 

AM 

mm 

ZVO-4  AU 

r- 

VO 

X 

-/ 

u» 

VO 

Uir-p-p- 

UI 

• 

»-4 

«n*rac3uj*x! 

X 

at 

4. 

acx*nrzaico 

at 

IO 

•■W 

«-i 

— « 

-4 

UI 

H-4Hp40H 

z 

OQOOOO 

1 1 1 1 1 

♦- 

r» 

UI 

a* 

vi  \ in  ru  co  0*4/1 

z 

■or»ai>rf"»*'6 

»■* 

A-l 

ui 

•-<ViU/*UU/»-4 

I- 

O' 

r»  cu-ooo*n* 

»— • 

• 

ONffltnoH 

< 

o 

X 

^urrunrupo*** 

■ 

CO 

r» 

— • 

a 

UJ 

p- 

O 

ui  ac»  < 

3 

•■u 

oxuu  ui 

UI 

VO 

a 

Zu*>o  at 

z 

< 

<Q  'ua 

o 

UI 

Ui 

*-  itzavi 

p- 

CO 

at 

«QiU<ZO 

a 

z 

m 

a. 

<>-OUi^<-> 

UI 

• 

API 

>1»>O*0LX 

p- 

ao 

<x 

VO 

-1 

*-* 

UI 

u- 

ac 

A*v 

H*u«4*^orw-* 

at 

^6 

vr 

oououoo 

O 

Ut 

rt* 

till  II 

P— 

— « 

immiKuumiUi 

o 

sro*L/*'*i"i<'iAn 

■T\r>{j'<jutrr-^r 

.*• 

mvoiww-u*-© 

VO 

o*Ao*u«-»r<i«u 

o* 

z 

MA«u«ur>iu»>rvo 

m* 

Ui 

• •••••• 

• • 

UI 

U1WAU4U4<-I«rt>4 

W 

z 

Ui 

zz 

UJ4« 

-J 

Z4WWUI 

o 

<* 

Z<-«4ZXZ 

CO 

at 

<J  Ui  HI  -a. 

■ ■■  • 

P- 

z-«au; 

z 

uu»"uvjua 

Ui 

uiui  at  ^ 

u 

Figure  lid. 


Histogram  and  sample  statistics  for  combined  waiting-time 
data  £ith  zero  waiting  times  removed)  in  the  simulated, 
correlated  MQxMD/l  queue;  t * 0.90,  g = 0.0,  p » 0.0, 
1/E(S)  * 4.0,  i.e.  all  Wn(j)'s  for  j = 1,...,500  and 
n = 21,000(1000)30,000. 
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Figure  lie.  Histogram  and  sample  statistic  for  combined  waiting-time 
data  (with  zero  waiting  times  removed)  in  the  simulated, 
correlated  MDcMD/1  queue;  t = 0.99,  3 * 0.0,  p = 0.0, 
1/E(S)  * 4.0,  i.e.  all  Wn(j)'s  for  j = 1,...,500  and 
n =■  21,000(1000)30,000. 
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CORRELATED  QUEUE  T = .99  R=.25  B=.50 


UNCGRRELATED  QUEUE  T = .99 


CORRELATED  QUEUE  T=.995  R=.25  B=.50 


VI.  VARIANCE  REDUCTION  IN  THE  SIMULATION 


A.  THE  MULTIPLE  CONTROL  VARIATES  METHOD 

Suppose  we  want  to  estimate  the  mean,  ux  of  a random 
variable  X and  suppose  V is  a zero-mean  random  variable, 
or  can  be  made  so  by  subtracting  from  V its  known  mean. 
Then 


(VI.  1) 


is  an  unbiased  estimator  for  ux  and 


(VI.  2) 


2 2 

where  ax  is  the  variance  of  X,  av  is  the  variance  of  V. 

2 2 2 

In  order  for  cr  to  be  a minimum  given  o ,a  and  cov(X,V], 
y xv 

the  optimal  choice  for  a is  readily  found  to  be 


(VI. 3) 


Cov[X,V] 


and  then,  substituting  back  in  (V.III),  the  minimum  variance 


2.2.  Cov  _ [X_i_Vj 

y * I* 

V 


2 2 
V1  - rx,v> 


where  r „ is  the  coefficient  of  linear  correlation  between 
x#  v 

X and  V.  Thus  the  variance  of  Y,  as  an  estimator  of  ux, 
is  always  less  them  that  of  X if  |cov(X,V)  | > 1.  In  a 
simulation  X would  be  the  usual  estimator  of  px,  and  Y would 
be  some  random  variable  generated  in  the  simulation  which  is 
correlated  with  x. 

2 

Generally , COV  (X,V)  andov  are  unknown,  and  we  can  use  in 
its  stead  an  estimate  derived  from  the  simulation.  This 
idea  of  reducing  the  variance  of  an  estimator  with  a control 
variable  is  readily  extended  to  a vector  V of  k control 
variates  V^,  . ..,  with  covariance  matrix  Q.  Now  let 

Y = X - a'V 

where  a'  » (a^,  a^)  is  a row  vector  of  coefficients, 

and  standard  results  from  multivariate  analysis  show  that 

- ox  + a'Qa  - 2a'R  . 

Here,  * cov[X,Vj].  The  optimal  choice  for  a is 

a*  ^ Q ^ 

(assuming,  of  course  that  Q is  non-singular,  for  otherwise  s^orne 

control  variates  would  be  redundant) . Then  the  minimal  value 
2 

of  0^  is  given  by 


Let  M be  the  joint  covariance  matrix  of  the  control  variates,  V 
and  X,  i.e. , 


R Q 


Then  it  cam  be  seen  that  with  the  optimal 


2 det  M 
y ~ det  Q 


For  further  details  see  Beja  (.1968)  and  Kleijnen  (19.74)  . 

B.  SELECTING  THE  CONTROL  VARIABLES 

In  order  to  simulate  the  mean  waiting  time  in  the  MDxMD/1 
queue,  several  control  variables  derived  from  the  known  pro- 
perties of  the  M/M/1  queue  and  the  service  and  interarrival 
processes  for  the  MD  x MD/1  queue  were  examined  and  combined 
to  be  the  multiple  control  variables.  Note  that  in  the  simula- 
tion the  M/M/1  queue  and  the  MD  x MD/1  queue  are  run  with  common 
exponential  variates,  so  that  the  service  and  interarrival 
times  in  the  M/M/1  queue  are  the  variables  from  which  the 
EARMA  type  service  and  interarrival  times  in  the  MD  x MD/1 
queue  are  generated.  It  is  well  known  that  in  a single-server 
FIFO  queue  the  difference  between  the  cumulative  interarrival 


1 


mmm 


times  and  cumulative  service  times  is  a good  control  for  Wn* 

The  additional  control  variables  were  meant  to  give  any 

greater  variance  reduction.  In  the  simulation  run  with 

p 0.25,  t * 0.50  and  8 * 0.50,  subroutine  CONVAR  was  used 

to  compute  the  correlation  between  waiting  time  Wn,  average 

waiting  time  Wn  in  the  MD  x MD/1  queue  and  several  quantities 

with  known  means  from  uncorrelated  queue,  i.e.,  the  M/M/1 

waiting  time,  the  M/M/1  average  waiting  time,  the  number  of 

arrivals  since  the  last  regeneration  point  (Wn  = 0) , sum 

of  service  time,  etc.  The  controls  for  and  W are  quite 

n n 

different  and  are  examined  separately. 

1.  Controls  for  the  Waiting  Time  Wr 

The  control  variables  from  the  uncorrelated  queue 

which  have  been  selected  as  being  the  most  correlated  to 

W are  the  number  of  arrivals  since  the  last  regeneration 
n 

(V^ , the  waiting  time  (V2) , the  average  of  the  previous  100 

interarrival  times  (V3) , the  average  of  the  previous  100 

interarrival  times  (V^) . The  correlation  matrix  for  these 

control  variables  with  W„  is  shown  in  table  5. 

n 

Since  the  mean  of  these  four  variables  are  not  zero, 
the  known  expected  value  must  be  subtracted  from  each  and 
the  estimator  is 


Y * X - o'  (V  - E (V) ) , 


where 


: 
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Note  that  even  though  V4,  the  sum  of  the  waiting  times 
in  the  MD  x MD/1  queue  is  an  average  of  correlated  random 


variables,  the  variables  are  exponential  and  their  average 


is  still  u . 

s 

A controlled  estimate  for  K'  was  obtained  using  V^, 

V2'V3/V4  an<i  resu^ts  are  shown  in  Table  6 and  Figure  14. 

The  500  replications  makes  it  possible  to  estimate,  as 

in  Table  5,  the  values  of  the  components  in  a*.  In  Figure 

14  the  second,  fourth,  ...  box  plots  are  for  the  sample  of 

500  controlled  W (j"’s,  j = 1,  ...,  500.  There  is  clearly 

n 

less  variability  than  in  the  box  plots  (first,  third,  etc.) 

for  the  uncontrolled  values  W (j),  j = 1,  ...,  500.  Successive 

n 

pairs  of  box  plots  are  for  different  values  of  n. 

Table  6 gives,  for  n = 61,000  (1,000)  70,000  statistics 

of  the  controlled  and  uncontrolled  W (j) , j = 1,  ...,  500 
• n 

samples.  The  values  S(MEAN)  are  to  be  compared.  Thus  for 
n = 70,000,  S (MEAN) , the  estimated  standard  deviation  of  Wn, 
is  0.06375,  which  is  to  be  compared  to  the  value  0.04632  for 
the  controlled  sample  average.  Thus  the  ratio  of  the  standard 
deviations  is  0.04632/0.06375  = 0.726,  and  the  variance 

reduction  is  (0.726)  = 0.528.  This  is  not  as  much  as  might 

' 

have  been  expected  from  the  multiple  controls,  but  it  reflects 

I 

the  difficulty  of  getting  further  control  variables  which  are 
not  themselves  highly  correlated  with  the  previous  control 
variables . 

2.  Controis  for  the  Sample-path  Average  WR 


Different  control  variables  are  needed  for  Wn  and 
Wn  since  the  former  includes  all  waiting  times  up  to  n and 


the  latter  is  local.  Thus  WR  will  only  be  correlated  with 

local  controls,  while  will  be  highly  correlated  with 

controls  which  contain  the  whole  history  of  the  sample  path. 

Thus  5 control  variables  from  the  M/M/1  queue  were  selected 

to  control  W time,  as  follows: 
n 

V. : number  of  zero  waiting  times  up  until  n in  the 
M/M/1  queue? 

V2:  cumulative  interarrival  time  in  the  M/M/1  queue; 
V3?  cumulative  service  time  in  the  correlated  queue; 
V4:  cumulative  servie  time  in  the  M/M/1  queue; 

V5:  mean  waiting  time,  Wn  from  the  M/M/1  queue. 

The  estimated  matrix  from  a simulation  run  with  500 
replications  at  p = 0.25,  6 » 0.5Cf,  t = 0.50  are  shown  in 
Table  7 . 

The  expected  values  of  V which  are  used  to  center  V 
are  as  follows: 


EfV^)  =»  n(l-t)  , where  n » no.  of  arrivals; 


n 

E(V2)  * £ EtX^)  , where  E(X.)  * mean  interarrival 

i»l  time  of  arrivals  in 

the  correlated  queue; 


e<v3) 


n 

£ E(SC. ) = njj^,  where  E(SC. ) * mean  service  time 
i»l  b of  arrivals  in  the 

MD  x MD/1  queue; 
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E(V4) 


e(v5) 


n 

l E(S.) 
i-1  1 


t 

T=t  us 


where  E(S.)  = mean  service 

time  of  arrivals 
in  the  M/M/1  queue; 


where  ps  = mean  service  time. 


The  results  of  the  simulation  run  with  500  replications 

and  p = 0.75,  t =*  0.75,  B * 0.50  are  shown  in  Table  8 and 

Fig.  15.  From  Table  6 we  see  that,  for  n * 70,000,  S (MEAN) 

for  Wn  with  and  without  control  is,  respectively  0.001128 

and  0.001709  and  we  recall  that  S(MEAN)  is  the  estimate  of 

the  standard  deviation  of  W . These  values  are  of  course  much 

n 

smaller  than  the  corresponding  values  for  Wn*  Further  the 

ratio  of  the  estimated  standard  deviations  is 

0.001128/0.001709  = 0.660035  and  the  variance  reduction  is 

2 

(0.660)  = 0.4356,  a rather  healthy  value.  The  variance 

reduction  shows  up  dramatically  in  Figure  15,  which  is  the 
analog  of  Figure  14.  It  is  also  clear  that  the  controlled. 

Wn's  are  symmetric,  possibly  normally  distributed,  and  that 
they  have  reached  steady  state. 


- • - 


Control 

Variable 


1.0 

0.1679 

0.4518 

-0.1195 

0.1450 

0.1679 

1.0 

0.6657 

-0.2611 

-0.3270 

0.4518 

0.6657 

1.0 

-0.2452 

0.0330 

-0.1195 

-0.2611 

-0.2452 

1.0 

0.3937 

0.1450 

-0.  327 

0.0330 

0.3937 

1.0 

Control  Variables 


Number  of  arrivals  since  last  regeneration  in 
M/M/1  queue; 

Waiting  time  of  M/M/1  queue; 

Average  of  100  previous  interarrival  times; 

Average  of  100  previous  service  times  in  MDxMD/1 
queue . 


Table  5.  Correlation  Matrix  of  Waiting  Ti~<ie  Wn  in  the 
MDxMD/1  Queue  and  the  Control  Variables; 
p » .25,  t * .50  and  8 3 0.50. 
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Table  6.  Estimated  Sample  Characteristics  of  Waiting  Time  Samples,  W (j),  j*i 
With  and  Without  Controls,  at  10  Different  Values  of  n.  n 
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Figure  14.  Box  plots  of  samples,  with  and  without  controls,  for  estimating  the  mean  waiting  time 

In  a correlated  queue  with  t = 0.75,  3 = 0.5,  p = 0.75.  The  first,  third nineteenth 

box  plots  are  of  Wni(j),  j=l,...,500;  for  ni=61,000,  n2=62,000,  ...  , m0*70,000, 

respectively.  The  second,  fourth twentieth  box  plots  are  of  the  Wni(j)  with  a 

multiple  control.  Note  that  the  controlled  Wn,(j)  may  be  negative.  1 


I 


VII.  CONCLUSION 


A.  In  simulation  the  queue,  the  box  plots,  of  the  sample 
of  waiting  times,  Wn(j),  and  average  waiting  time,  Wn(j), 
are  helpful  in  exploring  the  convergence  to  steady  state. 
Furthermore  the  box  plots  of  the  controlled  waiting  times, 
compared  with  the  uncontrolled  waiting  times , are  useful 
when  examining  just  how  much  variance  reduction  is  obtained. 
Note  that  the  variance  reduction  attained  for  WR  is  equivalent 
to  cutting  down  the  sample  size  by  at  least  a half. 

B.  An  approximation  for  the  mean  waiting  time  for  MDxMD/1 
queue  in  the  case  of  8 - 0.50  and  t -*■  1 is  found  to  be  equal 
to  (0.5  + ^lp5,P)  times  the  mean  waiting  time  of  M/M/1  queue. 
Thus  the  approximate  waiting  time  is 

*<»>  - I=t  X.  <°-5  + ' 

C.  The  multiple  control  variate  technique  is  useful  in  re- 
ducing  the  variance  of  the  waiting  time  and  mean  waiting  time 
when  simulating  the  correlation  queue.  This  is  done  by  simu- 
lating it  in  parallel  with  the  uncorrelated  queue  and  then 
using  the  data  from  the  uncorrelated  queue  to  control  the 
correlated  queue.  We  found  that  the  effectiveness  of  the 
variance  reduction  scheme  is  better  for  W , the  sample  path 

A 

mean  waiting  time,  them  for  the  waiting  time  W . 


first  r busy  periods,  so  that 


l X . = N 

i*l  1 


Now  select  a customer  at  random  from  the  first  N,  let  Ej 
be  the  event  that  the  customers  are  in  a busy  period  of  length 
j,  and  let  S be  the  number  of  customers  preceding  the  selected 
customer  in  his  busy  period.  Then 


for  k * 0,1, ... . , j-1 


18 


Since 


lira  N./r  * p. 
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